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Abstract 

We examine the first order Fermi acceleration on the presumption that supernova remnant shocks 
cross ambient magnetic fields with various angles. These oblique shocks accelerate particles more 
efficiently than the parallel shocks and elevate the maximum energies achievable by the particles. 
The primary cosmic ray spectrum is strongly dependent upon these energies. We also consider 
the dependence of the injection efficiency and of spectral indices on obliquity. When indices and 
absolute fluxes at 10 12 eV are given for six nuclear groups from balloon-borne data, each energy 
spectrum develops a smooth rigidity dependent knee structure. The resultant total spectrum also 
behaves similarly and fits well with ground-based experimental data up to several 10 17 eV. It is 
shown as well that the chemical composition changes significantly from lighter to heavier nuclei 
as the energies of particles exceed the knee region. Other predicted curves are compared with the 
experimental data which they reproduce rather well. 

PACS numbers: PACS 96. 40. De, 96. 50. Fm, 98. 70. Sa. 
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I. INTRODUCTION 



In recent years, many measurements of the primary spectrum and chemical composition 
of cosmic rays have been carried out. Some of them are performed directly by balloon- 
borne experiments up to about 100 TeV. Others are by indirect ground-based experiments 
ranging from sub- TeV to some 100 EeV. These experiments have mainly been done to study 
the knee region. Recent reports indicate that the all-particle spectrum of cosmic rays does 
not present a sudden steepening, but rather a gradual bent near the energy at the knee, 

£knee = (2 - 5) X 10 15 eV, from £-(2-6-2.8) to £-(3.0-3.2) g_ Thig value Qf ^ smaUer 

than about a factor of 2 compared to those of earlier reports M. The chemical composition 
of cosmic rays below -Eknee is directly measured and is almost constant in terms of the 
average logarithmic mass, < In A >. However, the data of < In A > above E^ nee fairly 
scatter, because they are obtained by indirect air shower observations. It is also noted 
that the value of < In A > varies depending on the hadronic interaction models assumed in 
simulations of extensive air showers 0. 

On the other hand, it is generally believed that cosmic rays in the energy range up 
to -Eknee are accelerated mostly by a supernova remnant (SNR) shock although there still 
exist some arguments that the origin of galactic cosmic rays (GCRs) is not yet settled [|J. 
When particles cross the shock front back and forth frequently, they obtain high energies 
proportional to the shock velocity, which is called the first order Fermi acceleration (for 
review, see ||). This stochastic acceleration and deformation of source spectra due to 
energy-dependent diffusion from the galaxy can explain the power-law spectrum up to the 
knee. This scenario predicts that the maximum energy E max , that particles can achieve by 
the shock acceleration mechanism, nearly equals i?knee- 

Actually E max is estimated to be ~ 3Z x 10 13 eV § at most, where Z is the atomic 
number, which is nearly two orders lower than -Eknee- However, this value is achieved by 
parallel shock acceleration where the normal of the shock front is parallel to the direction 
of the interstellar magnetic field. This is different from nonparallel shocks, where the angle 
«i between the shock normal and the field direction is not zero. These oblique shocks can 
accelerate particles much more efficiently |7|, §, For extremely oblique (quasiperpendic- 
ular) shocks, the values of -E max are increased by several orders of magnitude compared to 
those in parallel collisions. The details will be described in the next section. 
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In this paper we will explain the knee behavior of the energy spectrum with the accelera- 
tion mechanism as caused by oblique shocks. Our model is based on a simple idea as follows 
fT0|| . The ejected materials from a supernova explosion expand into the interstellar medium 
(ISM) driving a shock wave. While the expanding shell sweep up its own mass of the ISM, 
the shock acceleration may be active. The ejected mass (~ 10 M Q ) with large momentum 
moves out freely at constant velocity 10 8 — 10 9 cm/s during the free expansion phase. The 
shape of the shock front is supposed to be almost spherically symmetric. The directions of 
the interstellar magnetic field lines will be, over a wide range, nearly random rather than 
well aligned. It is therefore assumed that the field lines meet the shock front at random 
angles ati, the cosines of which, 77, are uniformly distributed. A more detailed discussion will 
be presented in section III. Since -E max is strongly dependent upon i], -E max values extend 
from E max (r] — 1) = E crit in parallel shocks to E max (r] = r] m m) = E cnt in quasiperpendicular 
ones. The uniform distribution of rj permits that particles with E, E crit < E < E cut , have 
less chances to reach high energies than those with E < E cr { t . 



There is an important question in acceleration at oblique shocks fTl] , |12| , |13| , however, 
namely that the rate of accelerated particles (injection efficiency) decreases drastically as 
the shock obliquity increases. On the other hand, the obliquity produces a harder energy 
spectrum |T3|, |T|, |TB|, |T7]]. Since particles with E < E CT - lt can be considered to be accelerated 
at shocks with the whole range of rj, we do not apply modifications to the given spectra. For 
E > -E cr it, both corrections of the injection efficiency and the change of the spectral index 
are taken into consideration, in addition to the chance probability of having 77. Their fluxes 
decrease and cause a gradual change of spectral indices. In terms of rigidity, -R cr it and -R cu t, 
corresponding to E cv - lt and -Ecut, are the same irrespective of Z. The flux of each elemental 
group starts to bend at different -E cr it. Then the total flux shows a slow change around the 
knee. 

We first choose the spectral indices and absolute fluxes for six nuclear groups at 1 TeV 
as an assumption so as to fit balloon-borne data. Next, based on our model, we calculate 
the energy spectrum up to 10 6 TeV for each group and obtain therefrom the total spectrum 
and energy dependence of the chemical composition directly. These results are compared 
with recent data of balloon and air shower experiments. 
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II. THE MAXIMUM ENERGY 



In this section we estimate the maximum energy of particles accelerated by oblique shocks 
at supernova remnants. If magnetic field lines are not parallel to the shock normal, an electric 
field is induced in the rest frame of the shock. Particles drift along the shock front gaining 
energies due to this field. In order to describe the interactions of particles with shocks 
concisely, the de Hoffmann- Teller frame in which the electric field vanishes globally is widely 
used For oblique shocks, reflection at the shock front is more important than diffusion 
in the downstream region (Kirk Sz Heavens 1989 [IS]), so that a more rapid acceleration 
can occur than for parallel shocks. The maximum energies of accelerated particles can be 
even a few orders of magnitude larger in the case of a quasiperpendicular shock compared 
to parallel ones. 

Suppose a plane shock wave propagating at a non-relativistic speed U\ with respect to 
the upstream medium. The shock normal intersects with the upstream magnetic field Bi at 
an angle et\. A schematic view of the shock in the upstream rest frame is shown in Fig. |l|. 
Parameters for the downstream region in its own rest frame, with magnetic field inclination 
/3 2 and shock velocity U 2 , can be expressed in terms of the upstream ones [pf| . We assume 
that the amplitude of the magnetic field perturbations is sufficiently small in the following 
calculation. 



A. The dependence of E max on magnetic field inclination 

There are some energy loss processes significant for the limitation of particle acceleration. 
For example, ionization, bremsstrahlung, and synchrotron radiation are conceivable. These 
processes for protons or other nuclei, however, are not so important within the energy range 
in which we are interested. The time scale of acceleration by oblique shocks at SNRs is 
regarded as sufficiently small compared to those of such loss processes. 

Unless the energy losses of particles are significant, the finite lifetime of the supernova 
blast wave as a strong shock mainly limits the maximum energy per particle. Effective shock 
acceleration can occur until the swept up mass reaches the ejected mass ~ 10 M & . The life- 
time of the shock t s h is estimated to be about thousand years for ejecta expanding at velocity 
~ 10 8 cm/s into a medium of average density 1 proton/cm 3 . As long as the characteristic 
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length for diffusion is much less than the radius of the shock, the plane approximation is 
valid. Balancing the acceleration time scale with t s h after Lagage & Cesarsky pi], we can 
write down the following expression for E max , 



E max = — — -U\eB x Z 

rex 



2 sin 2 ai r(cos 2 ai + -4 sin 2 ai 
cos at + — + — - — — 



(cos z «i + r z sin «ij 2 



where r is the compression ratio of the shock, i.e. r = U1/U2, and x is the square root of the 
ratio of km to /kj_, where the two k's are the parallel and normal components 

of the diffusion coefficient with respect to the magnetic field lines. 

Substituting the various constants with plausible values, we obtain E max as a function of 
r](= cosai) and x, 



£ max = 2.51 x 10 16 [eV] ( (—)(-^ 




30/iG/ V 3 P C 7 \ 107in / s , 
V 2 , r[ V 2 + ^(l-^)} 



-1 

■ ( 2 ) 

[i] 2 + r 2 (l — f] 2 )} 2 1 

The maximum energy for parallel shocks (77 = 1), E CT a can be estimated from Eq. (|2]) in 
the case of strong shocks (r = 4) with x = 30 and with typical values of £>i, R s h and Ui, 

E crit = 1.25Z x 10 14 [eV]. (3) 

Here we evaluate the numerical value 2.51 x 10 16 according to the recent papers. Namely, 
due to non-linear effects, the maximum energies are enhanced by the factor of 80/9 compared 



to ordinary values ||22|| . We will comment on the rather large value of B\ p3| , \L% in the next 
subsection. It is noted that the value 1.25 x 10 14 in Eq. (||) is larger than the usual ones 
~ 3 x 10 13 . 

It is also compared to Berezhko's values [P^| . When an enhancement factor of 80/9 is 
used, this value is larger by one order of magnitude, ~ 1 x 10 15 . Even if replacing the 
parameters of his Eq. (36), where the enhancement factor due to weak shock modification 
is 10/3 by our typical values mentioned above and after the normalized value of the field 3 
niG is corrected to 3 fiG |[25|| , we get ~ 1.4 x 10 15 . Thus in other words we may say that 
the value in Eq. (|3|) is only about ten times larger for field values used in a conventional 
derivation of E ma-X . Therefore numerical values in Eqs. (||) and (|]) are not inadequate even 
if x has a strong effect. 
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In order to see how E max (r]) varies with different x values, x = 10, 30, 100 are tested since 
the value of x is considered to lie between 10 and 100. Other parameters are the same as 
the above ones. Figure shows the r] dependence of E max for Z = 1: proton. From Eq. (0), 
the value of E m3uX is proportional to the particle charge Z. As shown in Fig. 0, the -E max 
value is strongly dependent on the field inclination rj. E max in quasiperpendicular shocks is 
larger than that in parallel shocks by two or three orders of magnitude for each value of x. 
In the case of x = 30, E cnt = 1.84 x 10 17 eV is about 1470 times the value E cr i t with Z = 1 
in Eq. (§. 



B. Magnetic field in the supernova remnants 

As we have mentioned above, we take a .Bi-value somewhat larger than the average 
strength of the interstellar magnetic field of magnitude of microgauss. It is widely known 
that SNRs are usually accompanied with magnetic fields. The sufficiently strong field is 
possibly present in the accelerated region around SNR-shocks because of rapid expansion of 
ejecta into the interstellar medium. There is some possibility [^] that cosmic ray streams 
amplify non-linearly magnetic fields up to 1 mG. Recent measurement also suggest that some 
galactic SNRs possess intense magnetic fields of the order from milligauss to sub-milligauss 



in the downstream region |27| t3 HI' H • These observed values are suppressed by a factor 
of 5 or so in the upstream, and hence our choice B\ = 30 /iG is quite relevant. 

There are some investigations related to the magnetic field geometry in SNRs. In young 
remnants, the ambient field lines are stretched and amplified to be predominantly parallel 
to the shock normal on account of Rayleigh- Taylor instabilities ||31||. On the other hand, 



for older SNRs perpendicular shocks likely originate from the compression of the ambient 



interstellar field. Observations of synchrotron radiation also support this morphology |32 



33| . The value of B\ of our choice and isotropic distribution of r] may be appropriate as far 
as we are concerned with the whole behavior of cosmic rays. 

Here we also comment on another important parameter, x, which describes field fluc- 
tuation. In papers of Jokipii (1987) [fjj and Ellison et al. (1995) |TJ|], they use £ = \/r g , 
which simply relates to x as x 2 = £ 2 + 1 (In their original papers, they use the notation 
r) instead of £ but here we use £ in order to distinguish it from our r](= cos«i)). £ ~ 1 
corresponds to the Bohm diffusion (strong scattering) limit and £ — > oo when k± <C K\\. 
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According to ref. fl3|| , the x value is limited by x < v/U±, where v is the particle velocity, 
and in quasi-perpendicular shocks x < tanai. Our choices x = 30 and rj m i n = 1/30, which 
will be adopted later, however, do not satisfy the latter inequality. The fixed x value is also 
used for any rj between rj min and 1. 



III. OUR MODEL 



Based on the result that the angular dependence of the maximum energy is quite signif- 
icant, we propose a simple model. The following three issues are examined: (1) the chance 
probability of 77 occurrence, (2) the dependence of the injection efficiency on rj and (3) the 
rj dependence of the spectral index. 



A. Chance probability of r/ 

The overall morphology of SNR-shocks is considered to be almost spherically symmetric 
even if the shape of each is somewhat distorted. The ambient disordered magnetic field lines 
over a wide range intersect the shock normal at various angles. Moreover, as we described 
in the section II. B, there exist SNRs with various ages whose magnetic fields cross the 
shock fronts with a wide variety of angles. Since both young and old SNRs are considered 
as sources of energetic cosmic rays, averaging over the various ages of SNRs and different 
shapes of their shock fronts, we assume a uniform distribution of rj in our model. 

The value of rj is restricted within the range between 77 min and 1 to satisfy the condition 
that the shock front should not move at the light velocity or above in the H-T frame. The 
probability in the width dt] is then 

f(rj)drj = ^ ( Vmin < rj < 1) , (4) 

-L ^min 

where 77 min = U\jc ~ 1/30 and (1 — r/mm) -1 is the normalization factor. The range where 
< rj < r] min is neglected because it occupies only 1/30 of the whole r\ region. Let £ be 
the energy of a particle with E CT - lt < £ < E cut . E CT - lt and E cnt are the maximum energies in 
parallel (rj = 1) and extremely oblique (quasiperpendicular: 77 = ?7 m i n ) shocks, respectively. 
A particular value of rj corresponds to an unique £ since E max (r]) is a monotonic decreasing 
function of rj as shown in Fig. 0. Only oblique shocks in the range from the f] to ?7 m i n can then 
serve to give £. The fluxes at £ should therefore be reduced by the factor (r]—r] min ) / (1— r] min ). 
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B. Injection efficiency 



Baring, Ellison and Jones JTT], Q developed test particle Monte-Carlo simulations with 
large angle scattering processes (abbreviated as LA). They showed that the efficiency for 
injecting thermal particles into an acceleration mechanism strongly depends on obliquity 
and varies inversely with at (@Bni hi their notation), in the case of high Mach numbers. 
They also checked the £ ( ~ our x) dependence of the injection efficiency, e(r/). When their 
results are applied to the present case x ~ 30, the injection efficiency relative to r\ = 1 is 
nearly equal to or less than 0.01 at r] — 0.6. Namely oblique shocks with rj from ?7 min to 0.6 
contribute very little to the acceleration of particles. 



Naito and Takahara [16], [17]] also attempted test particle simulations in oblique shocks. 
Instead of directly presenting e(r]), they showed the pitch angle (/x) distribution in both 



cases of LA and pitch-angle scattering (abbreviated as PA) in their figures 3 and 4 [17 



respectively, for rj = 1.0,0.75,0.5,0.25. Using the number of injected particles for these rj 
values, we can estimate e(rj) by integrating with respect to \i from —1 to 1 in the downstream, 
where the \x distribution of their simulation agrees well with the expected one The 
relative injection efficiency e(rj) changes little even if the number of accelerated particles in 
the upstream are added, because the numbers are less than those in the downstream and 
also the ratios to rj = 1.0 are similar to those in the upstream. 

For LA, e{r]) = 1, 0.73, 0.49, 0.26, whereas for PA e(rj) = 1, 0.81, 0.59, 0.31 are 
evaluated corresponding to T] = 1, 0.75, 0.5, 0.25, respectively. It is argued |35| that LA is 



adequate for ^ ~ 1, while PA is acceptable for ^ < 0.1. In our case, (^) 2 ~ \ = is 
somehow between LA and PA, and e(j]) is similar in both cases. Therefore, we may simply 
take, 

e(r?) = rj (5) 

The detailed reason why the above two works present very different results on e{rj), is not 
clear. The input for simulations, however, differs for the two approaches. In Ellison et al. 
fi~3| , |36| , injected particles are thermal, the energies of accelerated particles are up to about 
10 MeV, and U\ is slow, say 5 x 10 5 m/s. Ellison has also presented the application to the 
real SNR [37|j. He considered somewhat older objects which are in Sedov phase. On the 
other hand, in Naito and Takahara ||17|| , the Lorentz factors of injected particles are fixed 
at 2, the factors of accelerated particles are followed up to 10 5 , and the shock is relativistic: 
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U\ ~ 0.1c. In the light of the present treatments or parameters, such as JJ\ = 10 7 m/s and 
considering that the energies concerned are higher than 1 TeV, we assume Eq. (|^) here. 



C. Change of spectral indicies 

In this subsection, we consider the dependence of the spectral index on obliquity. It has 
been frequently pointed out that one of the advantages of the first order Fermi acceleration 
mechanism is that it gives a power-law energy spectrum. Namely one gets for the differential 
spectrum dJ/dE, if we take terms of fluid velocities up to first order into account, 

%« E - (™ th « = ^|)- <"> 

For a strong shock (r = 4), the differential spectral index is q = 2. This index is expected 
to be modified into the final index 7 which is observed. 7 is somewhat steeper than q by 
considering the energy-dependent leakage from the galaxy |j38| . We do not go into detailed 
discussions of this change in this paper, even though different modifications are used for 
different nuclear groups in the next subsection. 

This dependence of q values on rj is examined. This dependence is assumed to apply to 7 
in the same form. Although nonlinear back reaction effects of accelerated particles reduce q 
values by about 0.5 at higher energies even in parallel shocks ||39|| , only the index dependence 
on rj is considered here. The dependence has been studied in earlier works Here, we 

again refer to Naito and Takahara []TT5|, |T7[. In ref. [TTJ| they showed that q is practically 2 
at 77 = 1, and diminishes as 77 decreases except for 77 77 min being rather independent for the 
three cases ^ = 0.01, 0.03, 0.1. It is also shown @, for & = o.l in the case of LA, q is 2.08, 



slightly larger than 2, at 77 = 1, and almost linearly decreases as 77 decreases. Then q reaches 
a minimum of 1.5 at 77 ~ 0.15. For i] min < 77 < 0.15, q starts to increase rapidly. While for 
the PA case, q is nearly equal to 2.1 at 77 = 1. As 77 decreases, q decreases monotonically 
and reaches 1.0 as 77 approaches r/ m i n . In both cases the energy spectra become harder as 
shocks are more oblique. 

We choose the LA model except for 77 ps r] min and assume simply that q is proportional 
to 77 such as, 

q(r)) = ar] + b (1 < 77 < r7 min ), (7) 

where a = 0.68, b =1.41 and 77 min = This dependence of power indices on obliquity 
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enhances intensities of accelerated particles by oblique shocks. The present choice leads to 
an underestimate of this enhancement. For the energies of particles £ less than E CI it, the 
enhancement factor should be unity, because of q levelling off by summing up over the whole 
region of rj, i.e. from 1 to i] m [ n . 

For E CTit < £ < E cut , there is a one-to-one correspondence of r] to £ as mentioned above. 
Taking into account the injection efficiency given by Eq. (|5|), we introduce a correction 
factor 



£ 



-q(r)')+q 0?a) 



<(V) = T — 1 rf [ WT ) drf, (8) 

1— ?7min J VmiD \ Verity 

where r) a — (1 + r] min )/2 is the average value of rj. The factor —q(rf) + q(r) a ) leads to a 
deviation of indices between 1 ~ 77 min and r] ~ r] min and hence, b in Eq. (0) is cancelled. 
When £ = E CI n, that is to say rj — 1, £(77 = 1) = 1 due to the normalization factor 
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._ i. The 7] dependence of indices compensates partially the injection efficiency but not 

-L ??min 

completely. So Q{rj) is less than unity such as for r] min < r\ < 0.8, 0.7 < £(rf) < 0.8. 

In addition to C{rj), the chance probability of r] given by Eq. (S) should be combined. 
Since the probability is independent of ^"(77) , after integrating Eq. (^) we obtain a final 
correction factor: 

9(V) = C(V) (9) 

In Fig. 3, C(v), ?~ ??min and girj) are sketched as dashed, dotted and solid curves, respectively. 

-f ^min 

In summary, the differential energy spectra can be expressed as: 

= Cm g(v) (E ciit <E< E cut ) 

= (E cut < E), (10) 

where C and 7 are fixed according to six nuclear groups in the next section. E crit is given 
by Eq. (H) and E cnt = 1.84Z x 10 17 eV. From Eq. (flOf) , one can see that the curve of the 
flux is different iable changing from linear to a smooth bent above E CT a in a log-log diagram. 

The velocity of the shock front becomes superluminal in the H-T frame in the region 
< rj < r/min, which we have neglected. If particles are accelerated effectively by some other 
special mechanism, however, we should take this range into consideration. 
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IV. RESULTS 



We calculate differential ffuxes of various mass groups. The obtained fluxes are multiplied 
by E 2 ' 5 in order to clarify the change of the spectral indices. The multiplied flux for each 
chemical component is normalized to the experimental data at E = 1 TeV, where E is the 
total (mass plus kinetic) energy per particle, to give 

= <t>E-^ (11) 

£=lTeV 

where (3 = 7 — 2.5. Here we take as an assumption two different types of values of \I> and (3 as 
shown in Table |IJ. Parameters with subscript "HEGRA" have been taken from the analysis 
of the HEGRA air shower experiment HO] and slightly modified such that the heavy and 
very heavy nuclei (HVH) parts are separated into \1/ of sub-Fe and Fe with the ratio 0.125 
to 1.0 and with the same (3. The values with suffix "Ours" have been determined so as to 
fit well various experimental data around E = 1 TeV. 

We do not discuss in detail why some different values of \l/ and (3 are taken for groups of 
elements as in Table |. The most notable difference of these two models is the exponent (3 of 
the Fe spectrum. The resultant all-particle spectrum with "HEGRA" parameters shows a 
more rapid decrease for high energies. Since x ~ (B/5B) 2 which is estimated to be between 
10 and 100, for simplicity, we fix hereafter the value of x at about its geometrical mean of 
x = 30 in the calculation of Eq. ([!]). 



A. The energy spectra of primary cosmic ray components 

In Fig. |] we plot the calculated spectra of p and He against the kinetic energy E k per 
nucleon. The theoretical spectra with "Ours" parameters are expressed by solid curves 
and those with "HEGRA" parameters by dashed ones. The energy spectrum of p starts 
to bend due to the g factor in Eq. ([10|) at the energy: E kjCTit (p) ~ 1.25 x 10 5 GeV/n. 
This steepening of the spectral index is visible in the recent data of the Tibet air shower 
experiment [42[ ^3], ^4[. Here we cite the data of Tibet (2001) labelled proton dominant 



(PD) model in their paper. But their results from heavy dominant model are almost the 
same with those from PD. For He, E kjCTit (H.e) = Z(He)/A(He) x E crit (p) — m p ~ 6.3 x 10 4 



GeV/n. JACEE data [46] seem to be consistent with this £7fc iCr it(He), while RUNJOB's fluxes 



HU are somewhat lower but their highest E k does not reach _Efc iCrit (He). 
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We also present spectra of the CNO, middle, and Fe-group in Fig. |. "middle" corre- 
sponds to the NeMgSi-group and the Fe-group to our "sub-Fe" plus "Fe" in Table |. The 
two types of curves are defined by the same notation as in Fig. f|. The fluxes of middle 
and those of Fe-group have been downscaled by factors of 10 and 100 with respect to the 
original values, respectively. The slow changes of spectral indices of these groups above 
Ek^cnt — 6.3 x 10 4 GeV/n for CNO and 6.2 x 10 4 GeV/n for middle are not clearly seen from 
the data, because it is difficult to obtain such data by direct observations. 

For Fe, E kyCrit ~ 5.8 x 10 4 GeV/n, the JACEE point with the highest energy (~ 10 5 
GeV/n) shows a large value but with a large error. On the other hand, the recent Tibet 
data around ~ 10 6 GeV/n [[|7j] show that the spectrum is constant or has only a small 
bent in the diagram of Fig. |^. This favors Ours choice over HEGRA's. 

We plot the theoretical curves of spectra for various nuclear groups with Ours choice of 
parameters versus the total energy per particle in Fig. ^. One can see that the spectra of 
lighter nuclei start to bend at lower energies. This is just due to Eqs. (Q) and (|10|). Namely, 
the changes occur at a fixed rigidity of 1.25 x 10 5 GV as shown in Eq. (|3|). Summing the 
fluxes of all these groups, we obtain the flux of all particles which is plotted as "total" in 
this figure. Since the spectra of each group are smooth (differentiable) as mentioned above, 
the total spectrum should also be smooth. The knee behavior in our model comes from the 
^-dependence of E mSuX due to the oblique shock acceleration expressed in the second line of 
Eq. (|r0"p. This is different from the HEGRA analysis jKJ, Biermann J5"D[] , Stanev et al. |5]J 



and Erlykin and Wolfendale |52| . In the first paper, the knee is explained by introducing an 
artificial break of slope at a fixed rigidity. The authors of the next two papers, interpret the 
knee as the superposition of sources from different phases of SNRs. In the last reference, the 
authors argue that a sudden steepening at the knee and the explosion of a single, recent, 
nearby supernova is most likely to explain cosmic rays around the knee. 

In Fig. [7| the two curves of the total flux, i. e. with Ours (the same as in Fig. ^) and 
HEGRA choices of parameters, are compared with experimental data. As seen in this figure, 
Ours curve fits well Tibet [47j and Akeno data up to several times 10 8 GeV, while the 



curve with HEGRA parameters is in good agreement with DICE p3L RH] and CASA-MIA 



|56| , in which both E and \1/ are about 30% lower than the Akeno data. The KASCADE 
experiment presented a preliminary result with a four component assumption (H, He, C, 
and Fe) |)7| . It is shown that the position of the knee depends on rigidity. That is also 
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one of our conclusions. Anyway, the total flux can be explained at least up to 10 8 GeV by 
the present model. Above that energy, the solid curve falls below the AKENO data. Other 
sources than SNRs in Galaxy may be needed. 



using the least square method. For HEGRA and Ours choices 7 = 2.713 ± 0.006 and 
2.664 ±0.002, respectively, between 1.5 x 10 5 GeV and 1.5 x 10 6 GeV; and 7 = 3.248 ±0.012 
and 3.092 ± 0.008, respectively, between 6 x 10 6 and 6 x 10 8 GeV are obtained. These 
values depend on the adopted energy range; but they can be compared with experimental 
data. For example, Tibet 7 JTJ gives 2.60 ± 0.04 below and 3.00 ± 0.05 above -Eknee and 
CASA-BLANCA [§ 2.72 ± 0.02 and 2.95 ± 0.02. 

B. The energy dependence of chemical composition 

Using the results of energy fluxes, we examine how the chemical composition of cosmic 
ray particles changes with energy. In Fig. |] the flux ratio of p plus He to all particles is 
plotted against the energy per particle. Theoretical curves suggest a change with increasing 
energies from mixed composition to heavy particle dominance. The curve with HEGRA 
parameters lies almost within the uncertainties of experiments quoted here. The relative 
abundance of light elements, however, tends to increase again beyond the energy of several 
times 10 6 GeV from experiments of DICE and HEGRA. 

As a commonly used indicator of the composition, we take the average value of the 
logarithm of the mass number A, as 



where is the flux of species i which denotes each elemental group as shown in Table |. 
Our two curves of < In A > versus energy E are shown and compared with various data in 
Fig. [|. Both curves monotonically increase with energy, which is naturally confirmed by 
the present model. The HEGRA parametrization gives lower values than Ours in the whole 
energy range. 

In Fig. H experimental data scatter and have large uncertainties especially near and 
above the knee. A similar trend is seen in Fig. |9]. This is due to the fact that ground- 



We express the two curves in Fig. |7| in single power laws below and above -Eknce by 



<\nA> 




(12) 
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based measurements can study the composition only indirectly. The values of < In A > 
significantly depends on the experimental observables to be treated and on the hadronic 
interaction models to be used in simulation processes. Fairly different < In A > values are 
reported by various groups. Some authors [53| say that < In A > decreases in energy regions 



higher than the knee. Other authors || show that < In A > has even a dip near E^ nee . 
Most of the other data indicate that < In A > increases with E. Preliminary results of the 
KASCADE experiment |57| show a sharp increase in logarithmic mass above several times 



of 10 6 GeV. On the other hand, the data extracted from equi-intensity curves depicted by 
Chacaltaya II [B5f are somewhat closer to Ours. In the present stage we cannot reach a 



definite conclusion about whether our curves are in good agreement with data or not. 

V. CONCLUDING REMARKS 

We consider the acceleration caused by oblique shocks in which outer magnetic field lines 
cross the shock normal at any angle a.\. We calculate E m&K as a function of rj = cosai, 
where E max is the maximum energy that particles can achieve (see Eq. (^)). It is shown 
that -E max with extreme obliquity is more than 1400 times compared to that with parallel 
shocks (rj = 1). In the present model, rj is distributed uniformly and the rj dependence of 
the injection efficiency and the change of spectral indices are taken in to account. We choose 
the two parameter sets of power indices and absolute fluxes at 10 12 eV for various elemental 
groups shown in Table |. Then the spectra and compositions of cosmic rays are calculated 
up to much higher energies, say ~ 10 17 eV. It is shown that E cvlt = 1.25Z x 10 14 eV, where 
spectral curves start to bend over and continue up to E cnt ~ 1470 x E CT i t . These results fit 
the proton spectrum well and explain a smooth knee behavior around 3 x 10 14 eV (see Fig. 

As to the total flux, both curves in Fig. [7] exhibit a smooth knee behavior. The "HEGRA" 
choice fits well DICE and CASA-MIA data, while "Ours" choice reproduces well Tibet and 
Akeno data up to several times 10 17 eV. However, since our model is limited by -E max , for 
energies higher than ~ 10 17 eV, more energetic sources of cosmic rays should be necessary. 
Favorite candidates for the sources of cosmic rays far above the knee may be microquasers, 
pulsars, active galactic nuclei |66 , gamma-ray bursts (GRBs) |)7| or other hitherto unknown 



sources. 
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As for the composition, our value of < In A > increases with the energy per particle 
as shown in Fig. |9|. Data from various experiments scatter and have large uncertainties, 
particularly above -Eknee- Thus the situation remains uncertain. Recent results by Fly's 
Eye, Haverah Park, and HiRes groups provide some interesting data of composition in the 
high energy region above 10 8 GeV as summarized in ref. [ |BT] . Data of Haverah Park and 
HiRes show that light components are abundant towards 10 9 GeV, while a relatively heavy 
composition is suggested by Fly's Eye. These may also suggest some other origins apart 
from SNRs. 

In the present paper, we fixed x, Bi, R s h, r and U\ (i. e. ?7 m in) in Eq. (§) to evaluate 
the dependence of -E max on rj quantitatively and to obtain numerical values which could 
be directly compared with data. Generally, there are various shocks at SNRs where these 
quantities differ considerably depending on their ages, sizes and so on. However, as far as 
we are concerned with the whole or averaged behavior of cosmic rays which are observed, 
our choices and fixing of these quantities seem to be appropriate. We neglect the extremely 
quasi-perpendicular case in which the de-Hoffmann- Teller frame cannot be used. These are 
problems remaining to be solved. 
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FIG. 1: A schematic diagram of an oblique shock in the upstream rest frame where the magnetic 
field B\ intersects the shock normal at a nonzero angle a\. 
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FIG. 2: The maximum energy E max for a proton with three values of x versus magnetic field 
inclination rj. a\ is an angle between the magnetic field and the normal of the shock front. The 
E max at rj = 1 and rj = r] m i n , are the energies obtained by parallel shock and quasiperpendicular 
shock acceleration, respectively. For the nucleus with the atomic number Z, E mSLX should be 
multiplied by Z. 
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FIG. 3: Correction factors due to the injection efficiency plus the change of index power (dashed), 
chance probability (dotted) vs. rj. The final reduction factor is expressed by the solid curve. 
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FIG. 4: Comparison of calculated fluxes of p and He with experimental data. The horizontal axis 
represents kinetic energy per nucleon (GeV/n). The solid and dashed curves represent Ours and 
HEGRA choices of parameters, respectively as shown in Table |[ The original data are RUNJOB 



[41], Tibet etc. H2], Tibet [43], Tibet (2001) pL EAS-Top pL and for others see ref. [41]. 
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FIG. 5: Comparison of calculated fluxes of CNO-, middle, and Fe-groups with experimental data. 
The original fluxes of middle and of Fe-group are multiplied by 1/10 and 1/100, respectively. The 
solid and dashed curves represent the same notation as in Fig. 0. Experimental data are taken 



from RUNJOB's compilation [41], Tibet [47| and PDQ BACH pL 
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FIG. 6: The predicted fluxes of total, proton, He, and other nuclear groups in the case of "Ours : 
parameter choices shown in Table |. 
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FIG. 7: Comparison of all-particle spectrum versus total energy per particle between experimental 
data from various groups and the present calculation. The solid and dashed curves represent the 
same as in Fig. |. The data are Tibet g§, AKENO |, ||, DICE g |§, CASA-MIA @ and 
CASA-BLANCA Q (here only five data points with higher energies are cited in order to avoid 
overcrowding since other 16 points of lower energies mostly overlap with CASA-MIA data). Others 
are cited from RUNJOB's compilation [41|. 
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FIG. 8: Energy dependence of (p + a) /All. Two predicted curves which come from two choices of 
parameters are the same as in Fig. ^. Experimental data are as follows: HEGRA |5S[| , DICE [54], 



EAS-Top 5£], and others are Watson's compilation [60]. The large rectangle area shows direct 
measurement. 
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FIG. 9: Energy dependence of chemical composition of primary cosmic rays in term of < hiv4 >. 
The solid and dashed curves come from Ours and HEGRA choices of parameters, respectively 
as shown in Table § Experimental data are as follows: JACEE RUNJOB @, DICE §§ 



(X max +Muon), CASA-MIA g|, CASA-BLANCA (by QGSJET interaction model) Q, EAS-TOP 
KASCADE hadrons||, KASCADE (2001) (preliminary data) f57|, and Chacaltaya and Cha- 



caltaya II are cited from 64] and JBq ], respectively. About half of recent data from KASCADE 
(2001) and Chacaltaya II are plotted to avoid confusion. 
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TABLE I: Two choices of parameters for calculation at 1 TeV for the absolute value and index 
(3 in the form of Eq. (|ll|). ^ is in [GeV 1 - 5 m _2 s _1 sr -1 ]. 



Elements 


z 


A 




■ JTL JZjVj JX/\ 


^ uurb 




P 


1 n 


1 n 


^447 




ouuu 


97 


Hp 




4 n 


9087 


n 19 


99D0 


n 99 


CNO 


7.26 


14.5 


885 


0.17 


1030 


0.07 


Middle 


12.8 


25.8 


692 


0.11 


630 


0.07 


sub-Fe 


21.0 


45.0 


99 


0.11 


90 


0.07 


Fe 


26.0 


55.9 


790 


0.11 


720 


0.0 


total 






8000 




7670 
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